function ds = eleoms(~, s, mu, cexh)
alpha = atan2(-s(11), -s(13));
beta = atan2(s(12), -(s(11)*sin(alpha) + s(13)*cos(alpha)));

ds = zeros(14, 1);
alpha = atan2(-s(11), -s(13));
beta = atan2(s(12), -(s(11)*sin(alpha) + s(13)*cos(alpha)));
ds = zeros(14, 1);
ds(1) = s(4);
ds(2) = s(5)/s(1);
ds(3) = s(6)/(s(1)*sin(s(2)));
ds(4) = s(5)^2/s(1) + s(6)^2/s(1) - mu/s(1)^2 + s(7)*sin(alpha)*cos(beta);
ds(5) = -s(4)*s(5)/s(1) + s(6)^2/(s(1)*tan(s(2))) - s(7)*sin(beta);
ds(6) = -s(4)*s(6)/s(1) - s(5)*s(6)/(s(1)*tan(s(2))) + s(7)*cos(alpha)*cos(beta);
ds(7) = s(7)^2/cexh;
ds(8) = (s(9)*s(5) + s(11)*(s(5)^2 + s(6)^2) - s(12)*s(4)*s(5) - s(13)*s(4)*s(6))/s(1)^2 + ...
    (s(10)*s(6) + s(12)*s(6)^2*cos(s(2)) - s(13)*s(5)*s(6)*cos(s(2)))/(s(1)^2*sin(s(2))) - ...
    2*s(11)*mu/s(1)^3;
ds(9) = (s(10)*s(6)*cos(s(2)) + s(12)*s(6)^2 + s(13)*s(5)*s(6))/(s(1)*sin(s(2))^2);
ds(10) = 0;
ds(11) = -s(8) + s(12)*s(5)/s(1) + s(13)*s(6)/s(1);
ds(12) = (-s(9)  - 2*s(11)*s(5) + s(12)*s(4))/s(1) + s(13)*s(6)/(s(1)*tan(s(2)));
ds(13) = (-2*s(11)*s(6) + s(13)*s(4))/s(1) - (s(10) + 2*s(12)*s(6)*cos(s(2)) - ...
    s(13)*s(5)*cos(s(2)))/(s(1)*sin(s(2)));
ds(14) = -s(11)*sin(alpha)*cos(beta) + s(12)*sin(beta) - s(13)*cos(alpha)*cos(beta) - 2*s(14)*s(7)/cexh;
end